/*
    Copyright 2010-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Definitions of the loadart function that loads an ART dump file into the snapshot.
// $Id$

#include <iostream>
#include <string>
#include "preferences.h"
#include "biniostream.h"
#include "constants.h"
#include "snapshot.h"
#include "boost/lambda/lambda.hpp"

using namespace std;

template<typename T_stream>
T_stream read_entity(binifstream& f, T_stream) {
  T_stream temp;
  f >> temp;
  if(!f.good())
    throw T_stream();
  return temp;
}

/** Loads an ART dump file with the cells defined as particles. This
    shouldn't really be used, because it leads to funny artifacts. ART
    data should instead be exported with Yt. */
bool mcrx::Snapshot::loadart(const string& snapfn)
{
  using namespace boost::lambda;
 
  cout << "Loading ART dump files: " << snapfn << endl;
  
  time_ = 0.0;

  // Units in ART files are fixed, we convert them to the units we want immediately.
  lengthunit_=1;
  massunit_=1;
  timeunit_=1;

  G_=0;
  units_["length"] = "kpc";
  units_["mass"] = "Msun";
  units_["time"] = "yr";
  units_["temperature"] = "K";

  // The convert the filename stub to the actual file names for gas and stars.
  const size_t stubpos = snapfn.find("@");
  const string gasfn = string(snapfn).replace(stubpos, 1, "GZ");
  const string starfn = string(snapfn).replace(stubpos, 1, "SZ");

  // read gas data
  binifstream snap;
  snap.open(gasfn.c_str());
  if(!snap) {
    cerr << " Error opening ART gas file " << gasfn << endl;
    return false;
  }
  snap.bigendian();

  cerr << "Reading ART gas cells dump file " << gasfn << endl;

  typedef blitz::TinyVector<float, 3> vec3f;

  // The gas file consists of a number of fortran records containing
  // the particle data, so we can loop and add them one by one.
  // loop exits when the read function throws upon eof
  try {
    int i=0;
    while(true) {
      int dumi=read_entity(snap, int());
      // should be 44 bytes in this field
      assert(dumi==44);
      // get cell size in pc
      const float size=read_entity(snap, float())*1e-3;
      // pos is in kpc
      gas_particles.pos_.push_back(read_entity(snap, vec3f()));
      // velocity is in km/s -- need to convert later
      gas_particles.vel_.push_back(read_entity(snap, vec3f()));
      // density is in atoms/cm3
      const float rho=read_entity(snap, float());
      // push back mass which then is in units of H-atoms*kpc^3/cm^3 -- convert later
      gas_particles.m_.push_back(rho*size*size*size);
      // temp is in K
      gas_particles.temp_.push_back(read_entity(snap, float()));
      gas_particles.teff_.push_back(gas_particles.temp_.back());
      // metallicity is just the sum of the SNII and SNIa metallicities
      gas_particles.z_.push_back(read_entity(snap, float())+read_entity(snap, float()));
      
      // make up numbers for the fields we don't know about
      gas_particles.order_[i]=i;
      gas_particles.id_.push_back(i++);
      // the "smoothing length" is half the cell size
      gas_particles.h_.push_back(0.5*size);
      gas_particles.sfr_.push_back(0.0);
      gas_particles.pid_.push_back(gas_particles.id_.back());

      dumi=read_entity(snap, int());
      // should be 44 bytes in this field
      assert(dumi==44);
    }
  }
  catch(int&) {};
  snap.close();

  cout << "  Read " << gas_particles.n() << " gas cells\n";

  // read star data
  snap.open(starfn.c_str());
  if(!snap) {
    cerr << " Error opening ART star particles file " << starfn << endl;
    return false;
  }
  snap.bigendian();

  cerr << "Reading ART star particles dump file " << starfn << endl;

  // The star file consists of a number of fortran records containing
  // the particle data, so we can loop and add them one by one.
  // loop exits when the read function throws upon eof
  try {
    int i=0;
    while(true) {
      int dumi=read_entity(snap, int());
      // should be 72 bytes in this field
      assert(dumi==72);
      // get id
      nstar_particles.id_.push_back(read_entity(snap, int()));
      nstar_particles.order_[nstar_particles.id_.back()]=i++;
      // pos is doubles, in kpc
      nstar_particles.pos_.push_back(read_entity(snap, vec3d()));
      // velocity is doubles, in km/s -- need to convert later
      nstar_particles.vel_.push_back(read_entity(snap, vec3d()));
      // mass is in Msun (and is actually creation mass)
      nstar_particles.m_.push_back(read_entity(snap, double()));
      nstar_particles.mcr_.push_back(nstar_particles.m_.back());

      // files have age, so we push back tf as -age (in Gyr) and set snaptime to 0
      nstar_particles.tf_.push_back(-read_entity(snap, float())*1e9);
      // metallicity is just the sum of the SNII and SNIa metallicities
      nstar_particles.z_.push_back(read_entity(snap, float())+read_entity(snap, float()));

      // make up numbers for the fields we don't know about
      nstar_particles.pid_.push_back(gas_particles.id_.back());

      dumi=read_entity(snap, int());
      // should be 72 bytes in this field
      assert(dumi==72);
    }
  }
  catch(int&) {};

  cout << "  Read " << nstar_particles.n() << " star particles\n";

  // fix mass units of gas particles
  transform (gas_particles.m_.begin(), 
	     gas_particles.m_.end(), 
	     gas_particles.m_.begin(), 
	     _1 * units::convert("kg*kpc^3/cm^3","Msun") * constants::mp);

  // this sets the size of the star particles to the softening length
  // as defined in the preferences
  init_disk_bulge_halo();
  
  return true;
}

